# PEDL_3.R

# Contains: Code to generate Figure 3 and Table A3

# Load libraries
library(pacman)
pacman::p_load(ggplot2, dplyr)

# Set working directory (YOUR FILEPATH HERE)
setwd("")

#######################################################################
#######################################################################

# Import data
primary <- read.csv(file = "primary_cleaned.csv", header = TRUE)

#######################################################################
#######################################################################

# Clean arm variable
primary$arm.factor <- factor(primary$arm, levels = c("0", "125", "250", "375"))
primary$arm.dollar <- factor(ifelse(primary$arm == 0, "Control", 
                                            ifelse(primary$arm == 125, "$1.7",
                                                   ifelse(primary$arm == 250, "$3.4", "$5.1"))),
                                     levels = c("Control", "$1.7", "$3.4", "$5.1"))

# Create primary outcome variable
primary$intervention_weight <- primary$intervention_refills2*14.2

primary$lpg_wt_used <- primary$bline_lpg + 
  ((primary$intervention_refills2 - 1)*14.2) + (14.2 - primary$eline_lpg)

# Participants did not buy any LPG, and had a greater cylinder weight 
# at endline than baseline. These have been recoded to 0.
primary$lpg_wt_used <- ifelse(primary$lpg_wt_used < 0, 0, primary$lpg_wt_used)


# Table A3
# Calculate quantiles by grouping variable
q = c(.25, .5, .75)

primary %>%
  group_by(arm.dollar) %>%
  summarize(quant25 = quantile(intervention_weight, probs = q[1]), 
            quant50 = quantile(intervention_weight, probs = q[2]),
            quant75 = quantile(intervention_weight, probs = q[3]))



# Figure 3
f3 <- ggplot(primary, aes(x = intervention_weight)) +
  geom_bar(fill = "gray60") +
  geom_vline(data = ddply(primary, "arm.dollar", summarize, 
                          favg = median(intervention_weight, na.rm = TRUE)), 
             aes(xintercept=favg), color = "black", lty = 2) +
  facet_wrap(~arm.dollar, scales = "free") +
  xlim(c(-10, 120)) +
  scale_y_continuous(limits=c(0,50)) +
  xlab("Weight, kg") +
  ylab("Count") +
  theme_bw()

f3